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Equation of State in Relativistic Magnetohydrodynamics: 
variable versus constant adiabatic index 
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ABSTRACT 

The role of the equation of state for a perfectly conducting, relativistic magnetized 
fluid is the main subject of this work. The ideal constant T-law equation of state, 
commonly adopted in a wide range of astrophysical applications, is compared with a 
more realistic equation of state that better approximates the single-specie relativistic 
gas. The paper focus on three different topics. First, the influence of a more realis- 
tic equation of state on the propagation of fast magneto-sonic shocks is investigated. 
This calls into question the validity of the constant T-law equation of state in problems 
where the temperature of the gas substantially changes across hydromagnetic waves. 
Second, we present a new inversion scheme to recover primitive variables (such as 
rest-mass density and pressure) from conservative ones that allows for a general equa- 
tion of state and avoids catastrophic numerical cancellations in the non-relativistic 
and ultrarelativistic limits. Finally, selected numerical tests of astrophysical relevance 
(including magnetized accretion flows around Kerr black holes) are compared using 
different equations of state. Our main conclusion is that the choice of a realistic equa- 
tion of state can considerably bear upon the solution when transitions from cold to 
hot gas (or viceversa) are present. Under these circumstances, a polytropic equation 
of state can significantly endanger the solution. 

Key words: equation of state - relativity - hydrodynamics shock waves - methods: 
numerical - MHD 



1 INTRODUCTION 



Recent developments in numerical hydrodynamics have 
made a breach in the understanding of astrophysical phe- 
nomena commonly associated with relativistic magnetized 
plasmas. Existence of such flows has nowadays been largely 
witnessed by observations indicating superluminal motion in 
radio loud active galactic nuclei and galactic binary systems, 
as well as highly energetic events occurring in proximity of 
X-ray binaries and super-massive black holes. Strong evi- 
dence suggests that the two scenarios may be closely related 
and that the production of relativistic collimated jets results 
from magneto-centrifugal mechanisms taking pl ace in the in- 
ner r egions of rapidly spinning accretion disks l|Meier et al.l 
l200lf) . 

Due to the high degree of nonlinearity present in the 
equations of relativistic magnetohydrodynamics (RMHD 
henceforth), analytical models are often of limited appli- 
cability, relying on simplified assumptions of time inde- 
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pendence and/or spatial symmetries. For this reason, they 
are frequently superseded by numerical models that appeal 
to a consolidated theory based on finite difference meth- 
ods and Godunov-type schemes. The propagation of rel- 
ativistic supersonic jets without magnetic fiel d has been 
studie d, for instance, in the pione ering work of Ivan Puttenl 
1 19931'); Duncan fc Hughes! (|l994Tl and , subsequently, b y 

" 199Sl k lAlov et all i| 19991 ); 



Marti et all (|l997l);lHardee et al 



Mizuta et alj (2004) and references therein. Similar inves- 



tigations in presence of poloidal and toroidal magneti c 
fields have b e en carried on by iNishikawa et al.l (| 19971 ); 
Koidd dl997l); i Komissarovl (f l999) an d mor e recently by 
(|2005l ); lMignone et all (|2005l ). 
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The majority of analytical and numerical models, in- 
cluding the aforementioned studies, makes extensive use of 
the polytropic equation of state (EoS henceforth) , for which 
the specific heat ratio is constant and equal to 5/3 (for a 
cold gas) or to 4/3 (for a hot gas). H owever, the theory of 
relativistic perfect gases (|Svnge|[l957| ) teaches that, in the 
limit of negligible free path, the ratio of specific heats can- 
not be held constant if consistency with the kinetic theory 
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is to be requ ired. This was shown in an even earlier work by 
iTaubl (|1948T ). where a fundamental inequality relating spe- 
cific enthalpy and temperature was proved to hold. 

Although these results have been known for many 
decades, only f ew investigator s seem to have faced this im- 
portant aspect. iDuncan et alJ l|l996f ) suggested, in the con- 
text of extragalactic jets, the importance of self-consistently 
computing a variable adiabatic index rather than using a 
constant one. This may be advisable, for example, when 
the dynamics is regulated by multiple interactions of shock 
waves, leading to the forma tion of shock-heated regions in 
an initially cold gas. Lately, IScheck et all (|2002h addressed 
similar issues by investigating the long term evolution of 
jets with an arbitrary mixture of electrons, proton s and 
electron-positron pairs. Similarly, Mcli ani et al.l (|2004h con- 
sidered thermally accelerated outflows in proximity of com- 
pact objects by adopting a variable effective polytropic index 
to account for transitions from non-relativistic to relativis- 
tic temperatures. Similar considerations pertain to models 
of Gamma Ray Burst (GRB) engines including accretion 
discs, which have an EoS that must account for a combi- 
nation of protons, neutrons, electrons, positrons, and neu- 
trinos, etc. and must include the effects of electron degen- 
eracy, neutroni zation, photodisintegration, optical depth of 



neutrinos, etc. dPopham et al 



depth ot 

al .1120021 ; 



- .119991; iDi Matteo et ... _ _ 
iKohri fc Mineshigei l2002i ; iKohri et al l l2005h . However, for 
the disk that is mostly photodisintegrated and optically 
thin to neutrinos, a decent approximation of such EoS is 
a variable F-law with r = 5/3 when the temperature is be- 
low m e c 2 /kb and T — 4/3 when above m e c 2 /kb due to the 
production of positrons at high temperatures that gives a 
relativistic plasma (Broderick, McKinney, Kohri in prep.). 
Thus, the variable EoS considered here may be a reasonable 
approximation of GRB disks once photodisintegration has 
generated mostly free nuclei. 

The additional complexity introduced by more elabo- 
rate EoS comes at the price of extra computational cost since 
the EoS is frequently used in th e process of obtaining numer- 
ical solutions, see for example. iFalle fc Komissarovl |l996). 
Indeed, for the Synge gas, the correct EoS does not have a 
simple analytical expression and the thermodynamics of the 
fluid becomes entirely formulated in terms of the modified 

Bessel functions. 

Recently iMignone etafl (|2005al . MPB henceforth) in- 
troduced, in the context of relativistic non-magnetized flows, 
an approximate EoS that differs only by a few percent from 
the theoretical one . The advantage of this approximate EoS, 
earlier adopted by iMathewsl (|l97ll ). is its simple analytical 
representation. A slightly better approxi mation, ba s ed on an 
analytical expression, was presented bv lRvu et ail (|2006l ). 

In the present work we wish to discuss the role of the 
EoS in RMHD, with a particular emphasis to the one pro- 
posed by MPB, properly generalized to the context of rel- 
ativistic magnetized flows. Of course, it is still a matter of 
debate the extent to which equilibrium thermodynamic prin- 
ciples can be correctly prescribed when significant deviations 
from the single- fluid ideal approximation may hold (e.g., 
non-thermal particle distributions, gas composition, cosmic 
ray acceleration and losses, anisotropy, and so forth). Nev- 
ertheless, as the next step in a logical course of action, we 
will restrict our attention to a single aspect - namely the use 
of a constant polytropic versus a variable one - and we will 



ignore the influence of such non-ideal effects (albeit poten- 
tially important) on the EoS. 

In [J2] we present the relevant equations and discuss the 
properties of the new EoS versus the more restrictive con- 
stant F-law EoS. In Sj3] we consider the propagation of fast 
magneto-sonic shock waves and solve the jump conditions 
across the front using different EoS. As we shall see, this 
calls into question the validity of the constant T-law EoS 
in problems where the temperature of the gas substantially 
changes across hydromagnetic waves. In SJJ] we present nu- 
merical simulations of astrophysical relevance such as blast 
waves, axisymmetric jets, and magnetized accretion disks 
around Kerr black holes. A short survey of some existing 
models is conducted using different EoS's in order to deter- 
mine if significant interesting deviations arise. These results 
should be treated as a guide to some possible avenues of 
research rather than as the definitive result on any individ- 
ual topic. Results are summarized in i|5] In the Appendix, 
we present a description of the primitive variable inversion 
scheme. 



2 RELATIVISTIC MHD EQUATIONS 

In this section we present the equations of motion for rel- 
ativistic MHD, discuss the validity of the ideal gas EoS as 
applied to a perfect gas, and review an alternative EoS that 
properly models perfect gases in both the hot (relativistic) 
and cold (non-relativistic) regimes. 



2.1 Equations of Motion 

Our starting point are the relativistic MHD equations in 
conservative form: 



d_ 

dt 



( D \ 

m 
B 

\ E ) 



+ V 



/ Dv \ 

Wt^f 2 vv — bb + Ip t 
vB - Bv 



\ 



= 0, 



(1) 



/ 



together with the divergence-free constraint V B = 0, where 
v is the velocity, 7 is the Lorentz factor, w t = (ph+p + b 2 ) is 
the relativistic total (gas+ magnetic) enthalpy, pt = p + b 2 /2 
is the total (gas+magnetic) fluid pressure, B is the lab-frame 
field, and the field in the fluid frame is given by 



b a = 1 {v B,^+v l (v B)}, 

7 

with an energy density of 
\B? 



r 



+ (v-B) 



(2) 



(3) 



Units are chosen such that the speed of light is equal to one. 
Notice that the fluxes entering in the induction equation 
are the components of the electric field that, in the infinite 
conductivity approximation, become 

ft = —v x B . (4) 

The non-magnetic case is recovered by letting B — *• in the 
previous expressions. 



Equation of state in RMHD 3 




10 icr 



Figure 1. Equivalent T (top left), specific enthalpy (top right), 
sound speed (bottom left) and specific internal energy (bottom 
right) as functions of temperature © = p/p. Different lines corre- 
spond to the various EoS mentioned the text: the ideal T = 5/3- 
law (dotted line), ideal T = 4/3-law (dashed line), TM EoS (solid 
line). For clarity the Synge-gas (dashed-dotted line) has been 
plotted only in the top left panel, where the "unphysical region" 
marks the area where Taub's inequality is not fulfilled. 



The conservative variables are, respectively, the labora- 
tory density D, the three components of momentum rrik and 
magnetic field Bk and the total energy density E: 



D = 



PI 



= (Dhj + \B\ 2 )v k -(v-B)B k , 



(v-B) 



(•>) 
(6) 

(7) 



E = Dhj-p+L^- + LS±- 

The specific enthalpy h and internal energy e of the gas 
are related by 

h = 1 + e + ^ , (8) 
P 

and an additional equation of state relating two thermody- 
namical variables (e.g. p and e) must be specified for proper 
closure. This is the subject of the next section. 

Equations (0)-([7|) are routinely used in numerical codes 
to recover conservative variables from primitive ones (e.g., 
p, v, p and B). The inverse relations cannot be cast in 
closed form and require the soluti on of one or more non- 
linear equations. iNoble et al.l (|2006h review several methods 
of inversion for the constant T-law, for which pe = p/(T — 1). 
We present, in Appendix [X] the details of a new inversion 
procedure suitable for a more general EoS. 

2.2 Equation of State 

Proper closure to the conservation law (JT| is required in 
order to solve the equations. This is achieved by specifying 
an EoS relating thermodynamic quantities. The theory of 



relativistic perfect gases shows that the specific enthalpy is 
a functio n of the tem perature 9 — p/p alone and it takes 
the form (|Svngelll957l ) 



K 2 {l/Q) ' 



(9) 



where Ki and K3 are, respectively, the order 2 and 3 modi- 
fied Bessel functions of the second kind. Equation ([9} holds 
for a gas composed of material particles with the same mass 
and in the limit of small free path when compared to the 
sound wavelength. 

Direct use of Eq. © in numerical codes, however, 
results in time-consuming algorithms and alternative ap- 
proaches are usually sought. The most widely used and pop- 
ular one relies on the choice of the constant T-law EoS 



h= 1 + 



-e, 



(10) 



where F is the constant specific heat ratio. However, iTaubl 
( 1948) showed that consistency with the relativistic kinetic 
theory requires the specific enthalpy h to satisfy 



(h - G) (h - 48) > 1 . 



(11) 



known as Taub's fundamental inequality. Clearly the con- 
stant F-law EoS does not fulfill (lll|l for an arbitrary choice 
of r, while © certainly does. This is better understood in 
terms of an equivalent r eq , conveniently defined as 



h-1 
-1-0 



(12) 



and plotted in the top left panel of Fig. [T]for different EoS. 
In the limit of low and high temperatures, the physically 
admissible region is delimited, respectively, by r cq ^ 5/3 
(for 9^0) and r cq < 4/3 (for 9 -> 00). Indeed, Taub's 
inequality is always fulfilled when T ^ 4/3 while it cannot 
be satisfied for F ^ 5/3 for any positive value of the tem- 
perature. 

In a recent paper, iMignone et al.l l|2005al ) showed that 
if the equal sign is taken in Eq. (jTTJ , an equation with the 
correct limiting values may be derived. T he resulting Eo S 
(TM henceforth), previously introduced bv lMathewsl l|l97lf ). 
can be solved for the enthalpy, yielding 



or, using ph — p + pe + p in (|11[) with the equal sign, 



p 



lii (pe + 2p) e + 2 pe 



3 (pe + p) 



e + 1 3 



(13) 



(14) 



Direct evaluation of r cq using (|13[l shows that the TM EoS 
differs by less than 4% from the theoretical value given by 
the relativistic perfect gas EoS (|9]). The proposed EoS be- 
haves closely to the T = 4/3 law in the limit of high tem- 
peratures, whereas reduces to the V = 5/3 law in the cold 
gas limit. For intermediate temperatures, thermodynamical 
quantities (such as specific internal energy, enthalpy and 
sound speed) smoothly vary between the two limiting cases, 
as illustrated in Fig. [T] In this respect, Eq. (|13|l greatly im- 
proves over the constant F-law EoS and, at the same time, 
offers ease of implementation over Eq. @. Since thermody- 
namics is frequently invoked during the numerical solution 
of {TJ, it is expected that direct implementation of Eq. {13} 
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in numerical codes will result in faster and more efficient 
algorithms. 

Thermodynamical quantities such as sound speed and 
entropy are computed from the 2 nd law of thermodynamics, 



10 2 



dS ~ d\ogp , 



(15) 



(16) 



where S is the entropy. From the definition of the sound 
speed, 

2 _ dp 
ae s 

and using de = hdp (at constant S), one finds the useful 
expression 

re 

- - 9 h _ J h 

9 5/1-89 



r-law EoS , 



h ft. 



(17) 



3ft ft -6 



TM EoS . 



where we set ft = dh/dQ. In a similar way, direct integration 
of (|15l) yields S — k log a with 



with ft given by (|13j) . 



T-law EoS . 
TM EoS . 



(18) 



3 PROPAGATION OF FAST 
MAGNETO-SONIC SHOCKS 

Motivated by the previous results, we now investigate the 
role of the EoS on the propagation of magneto-sonic shock 
waves. To this end, we proceed by constructing a one- 
parameter family of shock waves with different velocities, 
traveling in the positive x direction. States ahead and be- 
hind the front are labeled with Uo and U\, respectively, and 
are related by the jump conditions 



v 3 [U] = [F(U)] 



(19) 



where v s is the shock speed and [q] — qi — qo is the jump 
across the wave for any qua ntity q. The set of j ump condi- 
tions l|19p may be reduced (|Lichnerowic j Il976l ) to the fol- 
lowing five positive-definite scalar invariants 



[J]=0, 

N] = o , 

[H] = 



'<_£_ 

J 2 



j2 + > + */2] 



[h/p] 



= 0, 



= 0, 



[ft 2 ] + J 2 



+ 2H [p] + 2 



= 0, 



where 



J = p-y-Jsivs - v x ) , 

is the mass flux across the shock, and 

T) = --(vB) + -^B x . 

P 7 



(20) 
(21) 

(22) 

(23) 
(24) 

(25) 
(26) 
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Figure 2. Compression ratio (top panels), internal energy (mid- 
dle panels) and downstream Mach number (bottom panels) as 
functions of the shock four- velocity "f s v s . The profiles give the so- 
lution to the shock equation for the non magnetic case. Plots on 
the left have zero tangential velocity ahead of the front, whereas 
plots on right are initialized with v y o = 0.99. Axis spacing is 
logarithmic. Solid, dashed and dotted lines correspond to the so- 
lutions obtained with the TM EoS and the F = 4/3 and T = 5/3 
laws, respectively. 



Here 7 S denotes the Lorentz factor of the shock. Fast or 
slow magneto-sonic shocks may be discriminated through 
the condition ceo > cti > (for the formers) or «i < a?o < 
(for the latters), where a — h/p — TL. 

We consider a pre-shock state characterized by a cold 
(po = 10 -4 ) gas with density p = 1. Without loss of gen- 
erality, we choose a frame of reference where the pre-shock 
velocity normal to the front vanishes, i.e., v x o = 0. Notice 
that, for a given shock speed, J 2 can be computed from the 
pre-shock state and thus one has to solve only Eqns. (|21|| - 
(123. 



3.1 Purely Hydrodynamical Shocks 

In the limit of vanishing magnetic field, only Eqns. ([23} and 
(|24|) need to be solved. Since J 2 is given, the problem sim- 
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plifies to the 2x2 nonlinear system of equations 
H _ n 



J 2 + 



[P] 



hi hp 
pi pa 



[h/ P ] 



(27) 
(28) 



We solve the previous equations starting from v a — 0.2, 
for which we were able to provide a sufficiently close guess to 
the downstream state. Once the pi and pi have been found, 
we repeat the process by slowly increasing the shock velocity 
v s and using the previously converged solution as the initial 
guess for the new value of v a . 

Fig. [2] shows the compression ratio, post-shock inter- 
nal energy e\ and Mach number Vi/c s i as functions of the 
shock four velocity i> s 7 s . For weakly relativistic shock speeds 
and vanishing tangential velocities (left panels) , density and 
pressure jumps approach the classical (i.e. non relativistic) 
strong shock limit at ^y a v a w 0.1, with the density ratio be- 
ing 4 or 7 depending on the value of F (5/3 or 4/3, respec- 
tively). The post-shock temperature keeps non-relativistic 
values (O <C 1) and the TM EoS behaves closely to the 
r = 5/3 case, as expected. 

With increasing shock velocity, the compression ratio 
does not saturate to a limiting value (as in the classical 
case) but keeps growing at approximately the same rate 
for the constant T-law EoS cases, and more rapidly for the 
TM EoS. This can be better understood by solving the 
jump conditions in a frame of reference moving with the 
shocked material and then transforming back to our origi- 
nal system. Since thermodynamics quantities are invariant 
one finds that, in the limit hi 3> ho m 1, the internal energy 
becomes ei = 71 — 1 and the compression ratio takes the 
asymptotic value 



Pi 



= 71 + 



71 + 1 

r - 1 



(29) 



when the ideal EoS is adopted. Since 71 can take arbitrarily 
large values, the downstream density keeps growing indefi- 
nitely. At the same time, internal energy behind the shock 
rises faster than the rest-mass energy, eventually leading to 
a thermodynamically relativistic configuration. In absence 
of tangential velocities (left panels in Fig. [2]), this transi- 
tion starts at moderately high shock velocities (7 s v s > 1) 
and culminates when the shocked gas heats up to relativis- 
tic temperatures (0 ~ 1 4- 10) for ^y a v a > 10. In this regime 
the TM EoS departs from the F = 5/3 case and merges on 
the r = 4/3 curve. For very large shock speeds, the Mach 
number tends to the asymptotic value (r— 1) _1//2 , regardless 
of the frame of reference. 

Inclusion of tangential velocities (right panels in Fig. [5J) 
leads to an increased mass flux ( J 2 oc 70) and, consequently, 
to higher post-shock pressure and density values. Still, since 
pressure grows faster than density, temperature in the post- 
shock flow strains to relativistic values even for slower shock 
velocities and the TM EoS tends to the F = 4/3 case at even 
smaller shock velocities (7 s w s > 2). 

Generally speaking, at a given shock velocity, density 
and pressure in the shocked gas attain higher values for lower 
F cq . Downstream temperature, on the other hand, follows 




the opposite trend being higher as V e 



5/3 and lower 



when F 



4/3. 



Figure 3. Compression ratio (top), downstream plasma /3 (mid- 
dle) and magnetic field strength (bottom) as function of the shock 
four- velocity -y a v a with vanishing tangential component of the 
velocity. The magnetic field makes an angle 7r/6 (left) and 7r/2 
(right) with the shock normal. The meaning of the different lines 
is the same as in Fig. [2] 



3.2 Magnetized Shocks 

In presence of magnetic fields, we solve the 3x3 nonlin- 
ear system given by Eqns. (|22[). ([23} and (|24|l . and di- 
rectly replace rji — rjoho/hi with the aid of Eq. (|21[) . 
The magnetic field introduces three additional parameters, 
namely, the thermal to magnetic pressure ratio (/3 = 2p/b 2 ) 
and the orientation of the magnetic field with respect to 
the shock front and to the tangential velocity. This is ex- 
pressed by the angles ct x and a y such that B x — \B\ cosa x , 
By — \B\ sin a x cos a y , B z — \B\sina x sina y . We restrict 
our attention to the case of a strongly magnetized pre-shock 
flow with (3 = 2p /bo = 10~ 2 . 

Fig. [3] shows the density, plasma j3 and magnetic pres- 
sure ratios versus shock velocity for a x = tt/6 (left panels) 
and a x — n/2 (perpendicular shock, right panels). Since 
there is no tangential velocity, the solution depends on one 
angle only (a x ) and the choice of a y is irrelevant. For small 
shock velocities ("y s v a < 0.4), the front is magnetically driven 
with density and pressure jumps attaining lower values than 
the non-magnetized counterpart. A similar be havior is found 
in classical MHD (| Jeffrey fc Taniuti I Density and 
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Figure 4. Density ratio (top), downstream plasma /3 (middle) 
and magnetic field strength (bottom) as function of "y s v s when 
the tangential component of the upstream velocity is vt = 0.99. 
The magnetic field and the shock normal form an angle 7r/6. The 
tangential components of magnetic field and velocity are aligned 
(left) and orthogonal (right). Different lines have the same mean- 
ing as in Fig. [2] 



magnetic compression ratios across the shock reach the clas- 
sical values around y s v s ~ 1 (rather than ^sV s ~ 0.1 as in 
the non-magnetic case) and increase afterwards. The mag- 
netic pressure ratio grows faster for the perpendicular shock, 
whereas internal energy and density show little dependence 
on the orientation angle a x . As expected, the TM EoS mim- 
ics the constant F = 5/3 case at small shock velocities. At 
"fsVs < 0.46, the plasma j3 exceeds unity and the shock starts 
to be pressure-dominated. In other words, thermal pressure 
eventually overwhelms the Lorentz force and the shock be- 
comes pressure-driven for velocities of the order of v s w 0.42. 
When fsVg > 1, the internal energy begins to become com- 
parable to the rest mass energy (c 2 ) and the behavior of the 
TM EoS detaches from the F = 5/3 curve and slowly joins 
the F = 4/3 case. The full transition happens in the limit of 
strongly relativistic shock speeds, ^y a Vs < 10. 

Inclusion of transverse velocities in the right state af- 
fects the solution in a way similar to the non-magnetic case. 
Relativistic effects play a role already at small velocities 
because of the increased inertia of the pre-shock state in- 
troduced by the upstream Lorentz factor. For a x — n/6 
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Figure 5. Density contrast (top), plasma /3 (middle) and mag- 
netic field strength (bottom) for vt = 0.99. The magnetic field is 
purely transverse and aligned with the tangential component of 
velocity on the left, while it is orthogonal on the right. Different 
lines have the same meaning as in Fig. f2] 



(Fig. [4| , the compression ratio does not drop to small val- 
ues and keeps growing becoming even larger (< 400) than 
the previous case when v t = 0. The same behavior is re- 
flected on the growth of magnetic pressure that, in addi- 
tion, shows more dependence on the relative orientation of 
the velocity and magnetic field projections in the plane of 
the front. When a y = n/2, indeed, magnetic pressure at- 
tains very large values (b 2 jb% < 10 4 , bottom right panel in 
Fig. [4| . Consequently, this is reflected in a decreased post- 
shock plasma /3. For the TM EoS, the post-shock properties 
of the flow begin to resemble the F = 4/3 behavior at lower 
shock velocities than before, 7 s u s ~ 2 3. Similar consid- 
erations may be done for the case of a perpendicular shock 
(a x = 7t/2, see Fig. [5]), although the plasma (5 saturates 
to larger values thus indicating larger post-shock pressures. 
Again, the maximum increase in magnetic pressure occurs 
when the velocity and magnetic field are perpendicular. 



4 NUMERICAL SIMULATIONS 

With the exception of very simple flow configurations, the 
solution of the RMHD fluid equations must be carried out 




Figure 6. Solution to the mildly relativistic blast wave (problem 
1) at t = 0.4. From left to right, the different profiles give den- 
sity, thermal pressure, total pressure (top panels), the three com- 
ponents of velocity (middle panel) and magnetic fields (bottom 
panels). Computations with the TM EoS and constant T = 5/3 
EoS are shown using solid and dotted lines, respectively. 



Figure 7. Solution to the strong relativistic blast wave (problem 
2) at t = 0.4. From left to right, the different profiles give den- 
sity, thermal pressure, total pressure (top panels), the three com- 
ponents of velocity (middle panel) and magnetic fields (bottom 
panels). Computations with the TM EoS and constant T = 5/3 
EoS are shown using solid and dotted lines, respectively. 



numerically. This allows an investigation of highly nonlinear 
regimes and complex interactions between multiple waves. 
We present some examples of astrophysical relevance, such 
as the propagation of one dimensional blast waves, the prop- 
agation of axisymmetric jets, and the evolution of magne- 
tized accretion disks around Kerr black holes. Our goal is to 
outline the qualitative effects of varying the EoS for some in- 
teresting astrophysical problems rather than giving detailed 
results on any individual topic. 

Direct numerical integration of Eq. fTjl has b een 
achieved using the P LUTO code ( Mignone et al.ll2007l ) in 
UJ\ JO]and HARM (|Gammie et al.ll2003r ) in 3131 The new 
primitive variable inversion scheme presented in Appendix 
[X]has been implemented in both codes and the results pre- 
sented in £|4.1I were used for code validation. The novel in- 
version scheme offers the advantage of being suitable for a 
more general EoS and avoiding catastrophic cancellation in 
the non-relativistic and ultrarelativistic limits. 



4.1 Relativistic Blast Waves 

A shock tube consists of a sharp discontinuity separat- 
ing two constant states. In what follows we will be con- 
sidering the one dimensional interval [0, 1] with a discon- 
tinuity placed at x = 0.5. For the first test problem, 
states to the left and to the right of the discontinuity are 
given by (p,p, B y , B z )l = (1,30,6,6) for the left state and 
{p,p,B y ,B z ) R = (1, 1,0.7,0.7) for the right state. This re- 
sults in a mildly relativistic configuration yielding a max- 
imum Lorentz factor of 1.3 ^ 7 ^ 1.4. The second test 



consists of a left state given by (p,p, B v ,B z )l = (1, 10 3 , 7, 7) 
and a right state (p,p, B y , B z )r = (1, 0.1, 0.7, 0.7). This con- 
figuration involves the propagation of a stronger blast wave 
yielding a more relativistic configuration (3^7^ 3.5). For 
both states, we use a base grid with 800 zones and 6 levels 
of refinement (equiv. resolution = 800 • 2 6 ) and evolve the 
solution up to t = 0.4. 

Computations carried with the ideal EoS with T = 5/3 
and the TM EoS are shown in Fig. [6] and Fig. [7] for the 
first and second shock tube, respectively. From left to right, 
the wave pattern is comprised of a fast and slow rarefac- 
tions, a contact discontinuity and a slow and a fast shocks. 
No rotational discontinuity is observed. Compared to the 
r = 5/3 case, one can see that the results obtained with 
the TM EoS show considerable differences. Indeed, waves 
propagate at rather smaller velocities and this is evident at 
the head and the tail points of the left-going magneto-sonic 
rarefaction waves. From a simple analogy with the hydrody- 
namic counterpart, in fact, we know that these points prop- 
agate increasingly faster with higher sound speed. Since the 
sound speed ratio of the TM and V = 5/3 is always less 
than one (see, for instance, the bottom left panel in Fig. (TJ, 
one may reasonably predict slower propagation speed for the 
Riemann fans when the TM EoS is used. Furthermore, this is 
confirmed by computations carried with T = 4/3 that shows 
even slower velocities. Similar conclusions can be drawn for 
the shock velocities. The reason is that the opening of the 
Riemann fan of the TM equation state is smaller than the 
F = 5/3 case, because the latter always over-estimates the 
sound speed. The higher density peak behind the slow shock 
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Figure 8. Jet velocity as a function of the Mach number for 
different values of the initial density contrast r\. The beam Lorentz 
factor is the same for all plots, -yj, = 10. Solid, dashed and dotted 
lines correspond to the solutions obtained with the TM EoS and 
the r = 4/3 and T = 5/3 laws, respectively. 



follows from the previous considerations and the conserva- 
tion of mass across the front. 



4.2 Propagation of Relativistic Jets 

Relativistic, pressure-matched jets are usually set up by 
injecting a supersonic cylindrical beam with radius r& 
into a uniform s tatic ambient medium (see, for instance, 
iMartf et al.lll997l ). The dynamical and morphological prop- 
erties of the jet and its interaction with the surrounding are 
most commonly investigated by adopting a three parameter 
set: the beam Lorentz factor 7t,, Mach number = vt/c s 
and the beam to ambient density ratio r\ = pb/p m - The 
presence of a constant poloidal magnetic field introduces a 
fourth parameter /3f, = 2pt/b 2 , which specifies the thermal 
to magnetic pressure ratio. 



4-2.1 One Dimensional Models 

The propagation of the jet itself takes place at the velocity 
Vj, defined as the speed of the working surface that sepa- 
rates shocked ambient fluid from the beam material. A one- 
dimensional estimate of Vj (for vanishing magnetic fields) 
can be derived from momentum flux ba lance in the frame of 
the working surface (jMarti et al.lll997h . This yields 



Vi 



lb 



T}h h /h r , 



1 + lb\Jrih b /h 

n 



(30) 



where and h m are the specific enthalpies of the beam and 
the ambient medium, respectively. For given 7^ and den- 
sity contrast 77, Eq. (|30p may be regarded as a function of 




Figure 9. Computed results for the non magnetized jet at t = 90 
for the ideal EoS (r = 5/3 and V = 4/3, top and middle panels) 
and the TM EoS (bottom panel), respectively. The lower and 
upper half of each panels shows the gray-scale map of density 
and internal energy in logarithmic scale. 
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Figure 10. Position of the working surface as a function of time 
for T = 5/3 (circles), T = 4/3 (stars) and the TM EoS (dia- 
monds). Solid, dotted and dashed lines gives the one-dimensional 
expectation. 
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Figure 11. Density and magnetic field for the magnetized jet at 
t = 80 (first and second panels from top) and at t = 126 (third 
and fourth panels). Computations were carried with 40 zones per 
beam radius with the TM EoS. 



the Mach number alone that uniquely specifies the pres- 
sure pb through the definitions of the sound speed, Eq. (117[) . 
For the constant T-law EoS the inversion is straightforward, 
whereas for the TM EoS one finds, using the substitution 
9 = 2/3 suffix, 



Pb = V 3 

where t„ 
tion 



t- 



(31) 

satisfies the negative branch of the quadratic equa- 



te 



15-6- 



.Ml 



+ t 24- 10- 



.Ml 



+ 9 = 0, 



(32) 



with t = tanhx. In Fig. [5] we show the jet velocity for in- 
creasing Mach numbers (or equivalently, decreasing sound 
speeds) and different density ratios r\ — 10~ 5 , 10 -3 , 10 _1 , 10. 
The Lorentz beam factor is 7b = 10. Prominent discrepan- 
cies between the selected EoS arise at low Mach numbers, 
where the relative variations of the jet speed between the 
constant F and the TM EoS's can be more than 50%. This 
regime corresponds to the case of a hot jet (O w 10 in the 
r\ = 10~ 3 case) propagating into a cold (O w 10 -3 ) medium, 
for which neither the F — 4/3 nor the F = 5/3 approxima- 
tion can properly characterize both fluids. 



4-2.2 Two Dimensional Models 

Of course, Eq. (|30|l is strictly valid for one-dimensional flows 
and the question remains as to whether similar conclusions 
can be drawn in more than one dimension. To this end we 
investigate, through numerical simulations, the propagation 
of relativistic jets in cylindrical axisymmetric coordinates 
(r,z). We consider two models corresponding to different 
sets of parameters and adopt the same computational do- 
main [0, 12] x [0, 50] (in units of jet radius) with the beam 



being injected at the inlet region (r ^ 1, z — 0). Jets are in 
pressure equilibrium with the environment. 

In the first model, the density ratio, beam Lorentz fac- 
tor and Mach number are given, respectively, by n = 10 -3 , 
jb = 10 and Mb = 1.77. Magnetic fields are absent. Inte- 
gration are carried at the resolution of 20 zones per beam 
radius using the relativistic Godunov scheme described in 
MPB. Computed results showing density and internal en- 
ergy maps at t = 90 are given in Fig. [5]for F — 5/3, V = 4/3 
and the TM EoS. The three different cases differ in several 
morphological aspects, the most prominent one being the 
position of the leading bow shock, z ~ 18 when V — 5/3, 
z w 48 for T = 4/3 and z r; 33 for the TM EoS. Smaller 
values of F lead to larger beam internal energies and there- 
fore to an increased momentum flux, in agreement with the 
one dimensional estimate (|30|l . This favors higher propaga- 
tion velocities and it is better quantified in Fig. [10] where 
the position of the working surface is plotted as a function 
of time and compared with the one dimensional estimate. 
For the cold jet (F = 5/3), the Mach shock exhibits a larger 
cross section and is located farther behind the bow shock 
when compared to the other two models. As a result, the 
jet velocity further decreases promoting the formation of a 
thicker cocoon. On the contrary, the hot jet (F = 4/3) prop- 
agates at the highest velocity and the cocoon has a more 
elongated shape. The beam propagates almost undisturbed 
and cross-shocks are weak. Close to is termination point, 
the beam widens and the jet slows down with hot shocked 
gas being pushed into the surrounding cocoon at a higher 
rate. Integration with the TM EoS reveals morphological 
and dynamical properties more similar to the F = 4/3 case, 
although the jet is ~ 40% slower. At t = 90 the beam does 
not seem to decelerate and its speed remains closer to the 
one-dimensional expectation. The cocoon develops a thin- 
ner structure with a more elongated conical shape and cross 
shocks form in the beam closer to the Mach disk. 

In the s econd case, we c ompar e models C2-pol-l and 
Bl-pol-1 of iLeismann et al.l (|2005l ) (corresponding to an 
ideal gas with T = 5/3 and F = 4/3, respectively) with 
the TM EoS adopting the same numerical scheme. For this 
model, 77 = 10 -2 , Vb = 0.99, Mb = 6 and the ambient 
medium is threaded by a constant vertical magnetic field, 
B z = \/2pb- Fig. [11] shows the results at t — 80 and 
t = 126, corresponding t o the final integration times shown 
in lLeismann et all (|21) < )5 ) for the selected values of Y. For the 
sake of conciseness, integration pertaining to the TM EoS 
only ar e shown and the reade r is reminded to the original 
work by ILeismann et ail (|2005l ) for a comprehensive descrip- 
tion. Compared to ideal EoS cases, the jet shown here pos- 
sesses morphological and dynamical properties intermediate 
between the hot (r = 4/3) and the cold (F = 5/3) cases. As 
expected, the jet propagates slower than in model Bl-pol-1 
(hot jet), but faster than the cold one (C2-pol-l). The head 
of the jet tends to form a hammer-like structure (although 
less prominent than the cold case) towards the end of the 
integration, i.e., for t > 100, but the cone remains more con- 
fined at previous times. Consistently with model C2-pol-l, 
the beam develops a series of weak cross shocks and outgo- 
ing waves triggered by the interaction of the flow with bent 
magnetic field lines. Although the magnetic field inhibits the 
formation of eddies, turbulent behavior is still observed in 
cocoon, where interior cavities with low magnetic fields are 
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Figure 12. Magnetized accretion flow around a Kerr black hole 
for the ideal T-law EoS with T = 4/3. Shows the logarithm of the 
rest-mass density in colour from high (red) to low (blue) values. 
The magnetic field has been overlayed. This model demonstrates 
more vigorous turbulence and a thicker corona that leads to a 
more confined magnetized jet near the poles. 



Figure 14. As in figure H"2l but for the TM EoS. This EoS leads 
to turbulence that is less vigorous than in the T = 4/3 model but 
more vigorous than in the V = 5/3 model. Qualitatively the TM 
EoS leads to an accretion disk that behaves somewhere between 
the behavior of the T = 4/3 and V = 5/3 models. 



4.3 Magnetized Accretion near Kerr Black Holes 
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Figure 13. As in figure [T2l but for T = 5/3. Compared to the 
r = 4/3 model, there is less vigorous turbulence and the corona 
is more sharply defined. 



formed. In this respect, the jet seems to share more features 
with the cold case. 



In this section we study time-dependent GRMHD numerical 
models of black hole accretion in order to determine the ef- 
fect of the EoS on the behavior of the accretion disk, corona, 
and jet . We study three mode l s simi lar to the models stud- 
ied by iMcKinnev fc Gammiej (|2004r ) for a Kerr black hole 
with a/M « 0.94 and a disk with a scale height (H) to ra- 
dius (R) ratio of H/R ~ 0.3. The constant F-law EoS with 
T = {4/3, 5/3} and the TM EoS are used. The initial torus 
solution is in hydrostatic equilibrium for the F-law EoS, but 
we use the F = 5/3 EoS as an initial condition for the TM 
EoS. Using the V = 4/3 EoS as an initial condition for the 
TM EoS did not affect the final quasi-stationary behavior 
of the flow. The simplest question to ask is which value of F 
will result in a solution most similar to the TM EoS model's 
solution. 

More advanced questions involve how the structure 
of the accretion flow depends on the EoS. The previ- 
ous results of this paper indicate that the corona above 
the disk seen in the simul ations l|De Villiers et al" I 120031 ; 
IMcKinnev fc Gammiej [2004) will be most sensitive to the 
EoS since this region can involve both non-relativistic and 
relativistic temperatures. The corona is directly involved 
is the production of a turbulent, m agnetized, thermal disk 
wind (Mc Kinney fc Naravan|[2006al fbl). so the disk wind is 
also expected to depend on the EoS. The disk inflow near 
the black hole has a magnetic pressure comparable to the gas 
pressure l|McKinnev fc Gammie]|2004h . so the EoS may play 
a role here and affect the flux of mass, energy, and angular- 
momentum int o the black hole. The magnetized jet asso- 
ciated with the Blandford & Ziiaiek ] solution seen in simu- 
lations (|McKinnev fc Gammiell2004 lMcKinnevll2006l ) is not 
expected to depend directly on the EoS, but may depend in- 



Equation of state in RMHD 11 



directly through the confining action of the corona. Finally, 
the type of field geometries observed in simulations that 
thread the disk and corona (Hiroseetal ] 120041 : iMcKinnevI 
120051 ) might depend on the EoS through the effect of the 
stiffness (larger F leads to harder EoSs) of the EoS on the 
turbulent diffusion of magnetic fields. 

Figs. If 21 If 31 and [14] show a snapshot of the accretion 
disk, corona, and jet at t ~ 1000GAf/c 3 . Overall the re- 
sults are quite comparable, as could be predicted since the 
r = {4/3, 5/3} models studied in iMcKinnev fc Gammia 
were quite similar. For all mo dels, the fiel d geom etries 
allowed are similar to that found in IMcKinnevI (|2005[ ). The 
accretion rate of mass, specific energy, and specific angular 
momentum are similar for all models, so the EoS appears to 
have only a small effect on the flow through the disk near 
the black hole. 

The most pronounced effect is that the soft EoS (r = 
4/3) model develops more vigorous turbulence due to the 
non-linear behavior of the magneto-rotational instability 
(MR1) than either the F = 5/3 or TM EoSs. This causes 
the coronae in the F = 4/3 model to be slightly thicker and 
to slightly more strongly confine the magnetized jet resulting 
in a slight decrease in the opening angle of the magnetized 
jet at large radii. Also, the F = 4/3 model develops a fast 
magnetized jet at slightly smaller radii than the other mod- 
els. An important consequence is that the jet opening angle 
at large radii might depend sensitively on the EoS of the ma- 
terial in the accretion disc corona. This should be studied in 
future work. 



5 CONCLUSIONS 

The role of the EoS in relativistic magnetohydrodynamics 
has been investigated both analytically and n umerically. The 
equatio n of state previously introduced by iMignone et al.l 
(2005a) (for non magnetized flows) has been extended to 
the case where magnetic fields are present. The proposed 
equation of state closely approximates the single-specie per- 
fect relativistic gas, but it offers a much simpler analyti- 
cal representation. In the limit of very large or very small 
temperatures, for instance, the equivalent specific heat ratio 
reduces, respectively, to the 4/3 or 5/3 limits. 

The propagation of fast magneto-sonic shock waves has 
been investigated by comparing the constant V laws to the 
new equation of state. Although for small shock veloci- 
ties the shock dynamics is well described by the cold gas 
limit, dynamical and thermodynamical quantities (such as 
the compression ratio, internal energy, magnetization and so 
forth) substantially change across the wave front at moder- 
ately or highly relativistic speeds. Eventually, for increasing 
shock velocities, flow quantities in the downstream region 
smoothly vary from the cold (r = 5/3) to the hot (F = 4/3) 
regimes. 

We numerically studied the effect of the EoS on shocks, 
blast waves, the propagation of relativistic jets, and magne- 
tized accretion flows around Kerr black holes. Our results 
should serve as a useful guide for future more specific stud- 
ies of each topic. For these numerical studies, we formu- 
lated the inversion from conservative quantities to primitive 
quantities that allows a general EoS and avoids catastrophic 
numerical cancellation in the non-relativistic and ultrarela- 



tivistic limits. The analytical and numerical models confirm 
the general result that large temperature gradients cannot 
be properly described by a polytropic EoS with constant 
specific heat ratio. Indeed, when compared to a more re- 
alistic EoS, for which the polytropic index is a function of 
the temperature, considerable dynamical differences arises. 
This has been repeatedly shown in presence of strong dis- 
continuities, such shocks, across which the internal energy 
can change by several order of magnitude. 

We also showed that the turbulent behavior of magne- 
tized accretion flows around Kerr black holes depends on the 
EoS. The F = 4/3 EoS leads to more vigorous turbulence 
than the r = 5/3 or TM EoSs. This affects the thickness of 
the corona that confines the magnetized jet. Any study of 
turbulence within the accretion disk, the subsequent genera- 
tion of heat in the coronae, and the opening and acceleration 
of the jet (especially at large radii where the cumulative dif- 
ferences due to the EoS in the disc are largest) should use 
an accurate EoS. The effect of the EoS on the jet opening 
angle and Lorentz factor at large radii is a topic of future 
study. 

The proposed equation state holds in the limit where 
effects due to radiation pressure, electron degeneracies and 
neutrino physics can be neglected. It also omits potentially 
crucial physical aspects related to kinetic processes (such 
as suprathermal particle distributions, cosmic rays), plasma 
composition, turbulence effects at the sub-grid levels, etc. 
These are very likely to alter the equation of state by ef- 
fectively changing the adiabatic index computed on merely 
thermodynamic arguments. Future efforts should properly 
address additional physical issues and consider more gen- 
eral equations of state. 
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APPENDIX A: PRIMITIVE VARIABLE 
INVERSION SCHEME 

We outline a new primitive variable inversion scheme that 
is used to convert the evolved conserved quantities into so- 
called primitive quantities that are necessary to obtain the 
fluxes used for the evolution. This scheme allows a gen- 
eral EoS by only requiring specification of thermodynamical 
quantities and it also avoids catastrophic cancellation in the 
non-relativistic and ultrarelativistic limits. Large Lorentz 
factors (up to 10 6 ) may not be uncommon in some astro- 
physical contexts (e.g. Gamma-Ray-Burst) and ordinary in- 
version methods can lead to severe numerical problems such 
as effectively div iding by zero and subtracti ve cancellation, 
see, for instance. iBernstein fc Hughes! (|2006f ). 

First, we note that the general relativistic conserva- 
tive quantities can be written more like special relativis- 
tic quantities by choosing a special frame in which to mea- 
sure all quantities. A useful frame is the zero angular mo- 
me ntum (ZAMO) obse rver in an axisymmetric space-time. 
See iNoble et al.l (|2006h for details. From their expressions, 
it is useful to note that catastrophic cancellations for non- 
relativistic velocities can be avoided by replacing 7 — 1 in 
any expression with (u a u a ) / '(7 + 1), where here u a is the 
relative 4-velocity in the ZAMO frame. From here on the 
expressions are in the ZAMO frame and appear similar to 
the same expressions in special relativity. 

Al Inversion Procedure 

Numerical integration of the conservation law ([TJ proceeds 
by evolving the conservative state vector U = (D, m, B, E) 
in time. Computation of the fluxes, however, requires veloc- 
ity and pressure to be recovered from U by inverting Eqns. 
([5])-([7]), a rather time consuming and c hallenging ta s k. For 
the constant-F law, a recent work by INoble et al.l (|2006h 
examines several methods of inversion. In this section we 
discuss how to modify the equations of motion, interme- 
diate calculations, and the inversion from conservative to 
primitive quantities so that the RMHD method 1) permits 
a general EoS; and 2) avoids catastrophic cancellations in 
the non-relativistic and ultrarelativistic limits. 

Our starting relations are the total energy density 



----- 

and the square modulus of Eq. ((6]l , 



(Al) 
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|m| 2 = (W + \B 



2\2, ,2 
M 



a2 

- J±_ (2W 

W 2 \ 



(A2) 



where S = m ■ B and W = Dh'y. Note that in order for 
this expression to be accurate in the non-relativistic limit, 
one should analytically cancel any appearance of E in this 
expression. Eq. (|A2[) can be inverted to express the square 
of the velocity in terms of the only unknown W: 

\B\ 2 ) + \m\ 2 W 2 



I ,2 S 2 (2W 
\v\ = - 



(W+\B\ 2 ) 2 W 2 
After inserting (|A3|) into (|A1[) one has: 



E= W-p + 



+ 



(A3) 



(A4) 



2 ' 2(\B\ 2 + W) 2 

In order to avoid numerical errors in the non-relativistic 
limit one must modify the equations of motion and sev- 
eral intermediate calculations. One solves the conservation 
equations with the mass density subtracted from the en- 
ergy by defining a new conserved quantity (E' = E — D) 
and similarly for the energy flux. In addition, operations 
based upon 7 can lead to catastrophic cancellations since 
the residual 7 — 1 is often requested and is dominant in the 
non-relativistic limit. A more natural quantity to consider 
is |i>| 2 or 7 2 |i;| 2 . Also, in the ultrarelativistic limit calcula- 
tions based upon 7(|i>| 2 ) have catastrophic cancellation er- 
rors when \v\ — *■ 1. This can be avoided by 1) using instead 
\u\ 2 = 7 2 |v| 2 and 2) introducing the quantities E' = E — D 
and W' — W — D, with W' properly rewritten as 

DW 2 



W 



+ X7 



(A5) 



1+7 

to avoid machine accuracy problems in the nonrelativistic 
limit, where \ = P e +P- Thus our relevant equations become: 

ISI 2 \B\ 2 \m\ 2 ~S 2 



E' = W' —p + 



\m\ 2 = (W + \B\ 



2(|B| 2 + W' + D) 2 ' 

: \u\ 2 

l + |u| 2 W 2 



(2W + \B[ 



(A6) 



(A7) 



where W = W + D. 

Equations (|A6f) and (|A7|) may be inverted to find W', 
p and \u\ 2 . A one dimensional inversion scheme is derived 
by regarding Eq. ()A6|) as a single nonlinear equation in the 



only unknown W' and using Eq. (|A7|I to express |it| 2 as a 
function of W' . Using Newton's iterative scheme as our root 
finder, one needs to compute the derivative 



dE 

aw 



1 



dp 
W 1 



(|S| 



s 2 



(A8) 



(\B\ 2 + W + Df ' 

The explicit form of dp/dW' depends on the particular EoS 
being used. While prior methods in principle allow for a 
general EoS, one has to re-derive many quantities that in- 
volve kinematical expressions. This can be avoided by split- 
ting the kinematical and thermodynamical quantities. This 
also allows one to write the expressions so that there is no 
catastrophic cancellations in the non-relativistic or ultrarel- 
ativistic limits. Assuming that p = p(x,p), we achieve this 
by applying the chain rule to the pressure derivative: 



( dp 
\dW> 



dp d\ dp dp 



+ 



(A9) 



d X dW dp dW ■ 

Partial derivatives involving purely thermodynamical quan- 
tities must now be supplied by the EoS routines. Derivatives 



with respect to W' , on the other hand, involve purely kine- 
matical terms and do not depend on the choice of the EoS. 
Relevant expressions needed in our computations are given 
in the Appendix. 

Once W' has been determined to some accuracy, the 
inversion process is completed by computing the velocities 
from an inversion of equation (JSJl to obtain 

W ( mfc + w Bk ) ■ (A10) 



Vk 



w + 



One then computes \ from an inversion of equation (| A5|) to 
obtain 



X = 



II" 



D\u\ 



(All) 



7 2 (l+7)7 2 ' 

from which p or pe can be obtained for any given EoS. The 
rest mass density is obtained from 

D 



7 



(A12) 



and the magnetic field is trivially inverted. 

In summary, we have formulated an inversion scheme 
that 1) allows a general EoS without re-deriving kinematical 
expressions; and 2) avoids catastrophic cancellation in the 
non-relativistic and ultrarelativistic limits. This inversion in- 
volves solving a single non-linear equation using, e.g., a one- 
dimensional Newton's method. A similar two-dimensional 
method can be easily written with the same properties, and 
such a method may be more robust in some cases since the 
one-dimensional version described here involves more com- 
plicated non-linear expressions. 

One can show analytically that the inversion is accurate 
in the ultrarelativistic limit as long as 7^e~achino f° r 7 an d 
P/(P7 2 )^ e machmc for pressure, where e m ac hine ~ 2.2 x 10~ 16 
for double precision. The method used bv lNoble et alj (|2006l ) 
requires 7;$e^achme/10 ^ ue t° the repeated use of the ex- 
pression 7 = 1/vl — v 2 in the inversion. Note that we 
use 7 = yjl + \u\ 2 that has no catastrophic cancellation. 
The fundamental limit on accuracy is due to evolving en- 
ergy and momentum separately such that the expression 
E — \m\ appears in the inversion. Only a method that 
evolves this quantity directly (e.g. for one-dimensional prob- 
lems one can evolve the energy with momentum subtracted) 
can reach higher Lorentz factors. A n example t est pr oblem 
is the ultrarelativistic Noh test in lAlov et akl (|l999 : ) with 
p = 7.633 x 1(T 6 , T = 4/3, 1 - v = 10" 11 (i.e. 7 = 223607) 
This test has p/(p7 2 ) ~ 1.6 x 1CP 16 , which is just below 
double precision and so the pressure is barely resolved in 
the pre-shock region. The post-shock region is insensitive to 
the pre-shock pressure and so is evolved accurately up to 
7 ~ 6 x 10 r . These facts are have been also confirmed nu- 
merically using this i nversion within HA RM. Using the same 
error measures as m lAlov et all (|l999h we can evolve their 
test problem with an even higher Lorentz factor of 7 = 10 7 
and obtain similar errors of <0.1%. 



A2 Kinematical and Thermodynamical 
Expressions 



The kinematical terms required in equation (|A9f) may be 
easily found from the definition of W , 



W' = Dhj -D = D(j-l) + X7 2 



(A13) 
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by straightforward differentiation. This yields 



dx 
dW 

and 

dp 



D 



,d\v\ 2 



d(l/7) _ I>yd|u| 5 



dW dW' 2 dW ' 

where 

2 [3W{W + \B\ 2 ) + \B\ A ] + \m\ 2 W : 



d\v\ 2 
dW 



W 3 



{W+\B 



2\3 



(A14) 



(A15) 



, (A16) 



is computed by differentiating (|A3|I with respect to W (note 
that d/dW' = d/dW). Equation (|A14|) does not depend on 
the knowledge of the EoS. 

Thermodynamical quantities such as dp/dx, on the 
other hand, do require the explicit form of the EoS. For 
the ideal gas EoS one simply has 



P(X,P) = — f—X. 



(A17) 



where \ = P e + P- By taking the partial derivatives of (|A17P 
with respect to x (keeping p constant) and p (keeping x 
constant) one has 



dp 
dx 



r - 1 



3T = °- 

dp 



(A18) 



For the TM EoS, one can more conveniently rewrite 
d) as 

MP + X-P) = (x~P){x + 2p-p), (A19) 

which, upon differentiation with respect to x (keeping p con- 
stant) yields 

|p^2 X + 2p-5 P 

Similarly, by taking the derivative with respect to p at con- 
stant x gives 

gg= 2x ~ 5p . (A21) 
dp 5p + 5 X -8p v ; 

In order to use the above expressions and avoid catas- 
trophic cancellation in the non-relativistic limit, one must 
solve for the gas pressure as functions of only p and x an d 
then write the pressure that explicitly avoids catastrophic 
cancellation as {XyP} ~ * 0. One obtains: 

2X(X + 2p) 



p(x,p) = 



(A22) 



5( X + p) + v/9 X 2 + ISpx + 25p 2 

Also, for setting the initial conditions it is useful to be able 
to convert from a given pressure to the internal energy by 
using 



P£(P,P) = \ [pH 3P = 

2 \ 2p+ y/9p 2 + 4p 2 



(A23) 



which also avoids catastrophic cancellation in the non- 
relativistic limit. 



'(fc+i) 



W 



where 



W 



i(k) 



f(W') 



df{W')/dW 



f{W') = W' — E' — p + 



2(\B\ 2 + W + DY 



(A24) 
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and df{W')/dW = dE'/dW is given by Eq. 
The iteration process terminates when the residual 
\W' (k+1) /W'(k) - l| falls below some specified tolerance. 

We remind the reader that, in order to start the iter- 
ation process given by (|A24|I . a suitable initial guess must 
be provided. We address this problem by initializing, at the 
beginning of the cycle, W' ' = W + — D, where W+ is the 
positive root of 



V(W, 1) = , 



(A26) 



and V(W, \v\) is the quadratic function 

V(W,\v\) = \m\ 2 -\v\ 2 W 2 +(2W+\B\ 2 )(2W+\B\ 2 -2E).(A27) 

This choice guarantees positivity of pressure, as it can be 
proven using the relation 

T(W,\v\) 



P = 



2(2W + \B\ 2 ) 



(A28) 



which follows upon eliminating the (S/W) 2 term in Eq. (| A2|l 
with the aid of Eq. (|A1|I . Seeing that V(W, \v\) is a con- 
vex quadratic function, the condition p > is equivalent to 
the requirement that the solution W must lie outside the 
interval [W_ , W+], where V(W±, \v\) = 0. However, since 
V(W,\v\) > V{W,1), it must follow that W+ > W+ and 
thus W+ lies outside the specified interval. We tacitly as- 
sume that the roots are always real, a condition that is al- 
ways met in practice. 



A3 Newton-Raphson Scheme 

Equation ()A6I) may be solved using a Newton-Raphson it- 
erative scheme, where the (k + l)-th approximation to the 
W' is computed as 



